function Mesh = process_mesh_2d(PDE, Nx, Ny, basis_type, mesh_type)
Mesh.basis_type = basis_type;
if 4 == nargin
    Mesh.mesh_type = [];
    Mesh.left = PDE.domain(1);
    Mesh.right = PDE.domain(2);
    Mesh.bottom = PDE.domain(3);
    Mesh.top = PDE.domain(4);
    Mesh.Nx = Nx;
    Mesh.Ny = Ny;
    Mesh.hx = (Mesh.right-Mesh.left)/Nx;
    Mesh.hy = (Mesh.top-Mesh.bottom)/Ny;
    [node, elem, Mesh.bde] = process_square_2d(Mesh.left, Mesh.right, Mesh.bottom, Mesh.top, Nx, Ny);
    Mesh.Np = (Nx+1)*(Ny+1);
else
    Mesh.mesh_type = mesh_type;
    switch mesh_type
        case "7shape"
            filename = '../mesh_7shape.msh';
        case "circle"
            filename = '../mesh_circle.msh';
        case "hole"
            filename = '../mesh_hole.msh';
        case "Lshape"
            filename = '../mesh_Lshape.msh';
        case "NACA"
            filename = '../mesh_NACA.msh';
        case "NACA0012"
            filename = '../mesh_NACA0012.msh';
        case "square_bd"
            filename = '../mesh_square_bd.msh';
        otherwise
            error("Invalid mesh type.");
    end
    fileID = fopen(filename, 'r');
    data = textscan(fileID, '%f %f %f %f', 'Delimiter', ' ', 'MultipleDelimsAsOne', true);
    fclose(fileID);
    Np = data{1}(1);
    Ne = data{2}(1);
    node = [data{1}(2:Np+1), data{2}(2:Np+1)];
    elem = [data{1}(Np+2:Np+Ne+1), data{2}(Np+2:Np+Ne+1), data{3}(Np+2:Np+Ne+1)];
    bdary = [data{1}(Np+Ne+2:end), data{2}(Np+Ne+2:end)];
    Mesh.bde = boundaryedges(elem, bdary);
    Mesh.Np = Np;
end
[Mesh.node, Mesh.elem, Mesh.bdn] = process_grid_interpolation(node, elem, Mesh.bde, basis_type);
Mesh.Ne = size(Mesh.elem,1);
Mesh.Nb = size(Mesh.node,1);
Mesh.Nlb = size(Mesh.elem,2);
end

%% process_square_2d
function [node, elem, bde] = process_square_2d(left, right, bottom, top, Nx, Ny)
[X, Y] = meshgrid(linspace(left, right, Nx+1), linspace(bottom, top, Ny+1));
node = [X(:), Y(:)];
Np = size(node,1);
k = (1:Np-Ny-1)';
cut = (Ny+1)*(1:Nx);
k(cut) = [];
elem = [k, k+Ny+1, k+1;
    k+1, k+Ny+1, k+Ny+2];
elem = sortrows(elem);
Ne = size(elem,1);
bd_bottom = [1:2*Ny:Ne;
    1:Ny+1:Np-Ny-1;
    Ny+2:Ny+1:Np]';
bd_right = [(Nx-1)*2*Ny+2:2:Ne;
    Np-Ny:Np-1;
    Np-Ny+1:Np]';
bd_top = [2*Ny:2*Ny:Ne;
    2*Ny+2:Ny+1:Np;
    Ny+1:Ny+1:Np-1]';
bd_top = sortrows(bd_top,1,"descend");
bd_left = [1:2:2*Ny-1;
    2:Ny+1;
    1:Ny]';
bd_left = sortrows(bd_left,1,"descend");
bde = [bd_bottom; bd_right; bd_top; bd_left];
end

%% boundaryedge
function bde = boundaryedges(elem, bdary)
% E_circle = [elem, elem(:,1)];
% E_bd = zeros(size(bdary,1),1);
% for i = 1:size(bdary,1)
%     for j = 1:size(elem,1)
%         if all(bdary(i,:) == E_circle(j,1:2)) || all(bdary(i,:) == E_circle(j,2:3)) || all(bdary(i,:) == E_circle(j,3:4))
%             E_bd(i) = j;
%         end
%     end
% end
% bde = [E_bd, bdary];
E_circle = [elem, elem(:,1)];
edges = [E_circle(:,1:2); E_circle(:,2:3); E_circle(:,3:4)];
[~, loc] = ismember(edges, bdary, 'rows');
E_bd = reshape(loc, size(elem,1), 3);
[rows, cols] = find(E_bd);
index0 = zeros(size(rows,1), 1);
for i = 1:size(rows,1)
    index0(E_bd(rows(i), cols(i))) = rows(i);
end
bde = [index0, bdary];
end

%% process_grid_interpolation
function [new_node, new_elem, new_bdn] = process_grid_interpolation(node, elem, bde, basis_type)
bdn = unique(bde(:,2:3), "stable");
switch basis_type
    case "P1"
        new_node = node;
        new_elem = elem;
        new_bdn = bdn;
    case "P1b"
        P = zeros(size(elem,1),2);
        for k = 1:size(elem,1)
            V = [node(elem(k,:),:)];
            P(k,:) = polygon_centroid(V);
        end
        new_col = (size(node,1)+1:size(node,1)+size(elem,1))';
        new_node = [node; P];
        new_elem = [elem, new_col];
        new_bdn = bdn;
    case "P2"
        v0 = reshape(elem', 3*size(elem,1), 1);
        v1 = reshape(([elem(:,2:3), elem(:,1)])', 3*size(elem,1), 1);
        totalEdge = [v0,v1];
        index1 = zeros(size(totalEdge,1),2);
        k = 1;
        for i = 1:size(totalEdge,1)
            if ~all(totalEdge(i,:) == [0, 0])
                for j = i+1:size(totalEdge,1)
                    if all(totalEdge(i,:) == flip(totalEdge(j,:)))
                        totalEdge(j,:) = [0, 0];
                        index1(k,:) = [j, i];
                        k = k+1;
                    end
                end
            end
        end
        uniqueEdge = totalEdge(~(totalEdge(:,1) == 0),:);
        P = (node(uniqueEdge(:,1),:) + node(uniqueEdge(:,2),:))/2;
        index0 = setdiff((1:size(totalEdge,1))', sort(index1(:)));
        index1 = sortrows(index1(~(index1(:,1) == 0),:));
        index2 = [find(totalEdge(:,1) ~= 0), (1:size(P,1))'];
        new_col = zeros(size(elem,1), 3);
        k = 1;
        for i = 1:3:size(totalEdge,1)-2
            new_col(k,:) = [i, i+1, i+2];
            k = k+1;
        end
        for k = 1:size(index1,1)
            idx = (new_col == index1(k,1));
            new_col(idx) = index1(k,2);
        end
        for k = 1:size(index2,1)
            idx = (new_col == index2(k,1));
            new_col(idx) = index2(k,2);
            idx = (index0 == index2(k,1));
            index0(idx) = index2(k,2);
        end
        [~,index3] = sortrows(bde(:,1));
        [~,index3] = sortrows(index3);
        index0 = index0(index3);
        index0 = index0 + size(node,1);
        new_node = [node; P];
        new_elem = [elem, new_col+size(node,1)];
        new_bdn = reshape([bdn'; index0'],[],1);
    otherwise
        error("Invalid basis type.");
end
end

%% polygon_centroid
function center = polygon_centroid(V)
n = size(V,1);
A = 0;
center_x = 0;
center_y = 0;
for i = 1:n
    j = mod(i,n) + 1;
    xi = V(i,1);
    yi = V(i,2);
    xj = V(j,1);
    yj = V(j,2);
    factor = xi * yj - xj * yi;
    A = A + factor;
    center_x = center_x + (xi + xj) * factor;
    center_y = center_y + (yi + yj) * factor;
end
center_x = center_x / (3 * A);
center_y = center_y / (3 * A);
center = [center_x, center_y];
end